Design, Synthesis, and Molecular Modeling Studies of a Novel Benzimidazole as an Aromatase Inhibitor

In this study, a series of novel 1,3,4-oxadiazole-benzimidazole derivatives were designed and synthesized. Their cytotoxic activities against five cancer cell lines, including A549, MCF-7, C6, HepG2, and HeLa, were evaluated by the MTT assay. The compounds 5b,c showed satisfactory potencies with much higher anticancer activity in comparison to the reference drug doxorubicin against the studied cancer cell lines. In vitro, enzymatic inhibition assays of aromatase (ARO) enzymes were performed. Molecular docking, molecular dynamics simulations, and binding free energy analyses were used to better understand the structure–activity connections and mechanism of action of the aromatase inhibitors. Two types of satisfactory 3D-QSAR (CoMFA and CoMSIA) models were generated, to predict the inhibitory activities of the novel inhibitors. Molecular docking studies were also carried out to find their binding sites and types of their interactions with the aromatase enzyme. Additionally, molecular dynamics simulations were performed to explore the most likely binding modes of compounds 5b,c with CYP19A1.


INTRODUCTION
Cancer is a multifactorial complex disorder, characterized by the uncontrolled growth and spread of abnormal cells. 1 On the basis of a WHO survey, cancer has been responsible for almost 9.6 million deaths worldwide in 2018 and there will be around 13.1 million estimated new cases by 2030. 2,3 After cardiovascular disease, cancer is the second leading cause of death. 4 Although there are many drugs for the treatment of cancer, the need for continued research for new anticancer molecules remains due to selectivity, potency, and drug resistance. 5 Breast tissue has an excess of the aromatase enzyme. 6 Aromatase inhibitors (AIs) are a well-known endocrine therapy method that works by suppressing estrogen production, which is important for aromatase activity. 7 AIs are divided into four generations based on their clinical development order: first, second, third, and fourth generation AIs. AIs are divided into two types: type I (steroidal) and type II (nonsteroidal) on the basis of their mechanism of action. 8 The most effective, selective, and least toxic aromatase inhibitors are the third-generation aromatase inhibitors (anastrozole and letrozole). In actual practice, third-generation inhibitors provide a better therapeutic benefit with practically full specificity. Letrozole and anastrozole are nonsteroidal derivatives that have a triazole ring which interacts with the heme prosthetic group of aromatase and competes with androgen substrates by inhibiting them competitively. 9 A benzimidazole ring, which is an important core in medicinal chemistry, constitutes structures found in current drugs. The ring structure is also an isostere of DNA bases that carry purine and pyrimidine cores, and it can as well be a purine antimetabolite. 5 It has been proven by various studies that benzimidazole derivatives with anticancer activity show their activities by different mechanisms through DNA intercalation, topoisomerase I and II inhibition, an androgen receptor antagonistic effect, PARP-poly inhibition, dihydrofolate reductase (DHFR), aromatase inhibition, and microtubule inhibition. 10−15 Various classes of therapeutically active benzimidazole derivatives are available in the market such as bendamustine, veliparp, carbendazim, and nocodazole.
In our previous studies, we synthesized a series of benzimidazole-oxadiazole derivatives which showed satisfactory anticancer activity with good IC 50 values against various cell lines. 16,17 As part of a continuing effort to find out and screen more active benzimidazole-oxadiazole derivatives, in the present study, a series of novel 4-phenyl-and 4-piperazinylsubstituted 2-{[5-(2-phenyl-1H-benzo[d]imidazol-6-yl)-1,3,4oxadiazol-2-yl]thio}-1-piperazin-1-ylethan-1-one derivatives were synthesized which were later confirmed by 1 H NMR, 13 C NMR, and mass spectral techniques. Then, their cytotoxic activities against various cancer cell lines including A549 (human lung adenocarcinoma epithelial cell line), HepG2 (human colon cancer cell line), MCF-7 (breast adenocarcinoma), C6 (breast adenocarcinoma), and HeLa (ovarian cancer cell line) were evaluated. Using the crystal structure of aromatase retrieved from the Protein Data Bank (PDB ID 3EQM) and Autodock Vina software, docking experiments corroborated the results obtained from biological screening by predicting probable binding interactions of the target compounds with aromatase active sites. Molecular dynamic simulations were applied to evaluate the stability in the aromatase active site to justify their inhibitory activity.

RESULTS AND DISCUSSION
2.1. Chemistry. Scheme 1 depicts the synthesis of the target compounds 5a−i. In the first stage, a 4-substituted benzaldehyde was heated with 3,4-diaminobenzoic acid in DMF and sodium metabisulfite to obtain the derivatives 1a−c. By using a Fischer esterification procedure, compounds 1a−c were transformed into the methyl ester compounds 2a−c. To prepare compounds 3a−c, appropriate solutions of compounds 2a−c in ethanol (95%) were treated with hydrazine hydrate. Compounds 4a−c were formed by reacting hydrazide derivatives 3a−c with carbon disulfide in ethanolic potassium hydroxide. To make the target compounds 5a−i, compounds 4a−c were reacted with acetylated piperazine derivatives in acetone in the final reaction step. Using current analytical techniques such as 1 H NMR, 13 C NMR, and HRMS, the structures of synthesized compounds 5a−i were studied.
The S−CH 2 (methylene) protons were detected as a singlet between 4.53 and 4.63 ppm, according to a 1 H NMR spectral analysis of compounds 5a−i. At 3.04−3.86 ppm, the protons of the piperazine moiety appear as a multiplet. The protons of the methoxy substituent provided a singlet signal at 3.84−3.86 ppm in the 1 H NMR spectra of compounds 5d−f with a 4methoxyphenyl group in the second position of the benzimidazole ring. The −OCH 2 protons were found at 4.02−4.11 ppm in the −OC 2 H 5 group of compounds 5g−I on the phenyl ring, and CH 3 protons were found at 1.24−1.35 ppm. Due to the H-5 proton, the benzimidazole proton was observed as a doublet of doublets in 1 H NMR spectra at around 7.68−7.87. The carbon atoms of the molecule have chemical shift values that are similar to those reported in the literature in 13 C NMR. The masses of the target compounds were validated by an HRMS analysis, which confirmed the estimated values.
2.2. Cytotoxicity Assay. In this study, all compounds were tested for antiproliferative activity using an MTT assay with  Table 1 summarizes the results of four independent studies as the mean IC 50 value (half-maximum inhibitory concentration). The benzimidazole derivatives 5b,c possessing 4-chlorobenzyl and 4-fluorobenzyl groups were the most potent compounds against the MCF-7 (IC 50 = 5.132 ± 0.211, 6.554 ± 0.287 μM) cell line, respectively. Furthermore, compounds 5a,I exhibited moderate activity against MCF7. Compound 5c was the most active agent against the HeLa cell line (IC 50 = 7.316 ± 0.276 μM). Furthermore, compound 5i exhibited the best activity against A549 and C6 cell lines. In addition, among the compounds 5a−i, compounds 5a,i exhibited anticancer activity similar to that of doxorubicin against MCF7 cell line with IC 50 values of 9.963 ± 0.312 and 9.056 ± 0.311 μM, respectively.
According to the findings, the 4-hydroxyphenyl group boosted anticancer activity in MCF-7, C6, and HepG-2 cancer cell lines more than the 4-methoxyphenyl and 4-ethoxyphenyl groups. The presence of chlorobenzyl and fluorobenzyl groups at the fourth position of the piperazine scaffold also increased anticancer activity, in comparison to 4-chlorophenyl and benzyl substituents.
2.3. Aromatase Inhibition Assay. An in vitro fluorescence-based test (Aromatase (CYP19A) Inhibitor Screening Kit (Fluorimetric), BioVision) was used to determine the inhibitory activity of the produced compounds against human aromatase. The target drug was dissolved in acetonitrile, and the results were compared to those for the reference compound, letrozole. According to the result, compounds 5b,c effectively caused 50% aromatase enzyme inhibition at 1.475 ± 0.062 and 2.512 ± 0.124 μM, respectively ( Table 2).
2.4. Molecular Docking Analysis. A molecular docking analysis of 5b,c, two of the most active compounds designed and synthesized in this study, was performed. For the molecular docking study, the 3D crystal structure of the aromatase enzyme was downloaded from the Protein Data Bank (PDB ID 3EQM, resolution 2.90 Å). Heteroatoms in the structure, except for the HEM structure, were removed. The active site amino acids of the protein were determined on the basis of the cocrystal 4-androstene-3,17-dione (ASD) natural ligand. Molecular docking validation was achieved by removing the ASD ligand found in the 3EQM crystal structure and redocking the same area. Between the cocrystal ligand and the binding mode obtained by docking, the RMSD was measured to be 0.46 Å. Ligand structures were minimized using a universal force field (UFF). As shown in Figure 1, compounds 5b,c gave almost the same geometry binding pose at the CYP19A1 active site. Compound 5b formed a conventional H bond and a Hem600 hydrophobic π-sulfur interaction between the NH of oxazole and Asp309 at the CYP19A1 active site. Other interactions of compound 5b are given in Table 3. Compound 5c, on the other hand, formed Glu493 with phenol −OH, conventional H bonds with oxazole NH and Asp309, and a hydrophobic π−sulfur interaction with Hem600. Electrostatic and hydrophobic interactions between CYP19A1 and 5c are detailed in Table 3. Also, the interactions and binding poses of compounds 5b,c and the reference compound letrozole used in the experimental study at the CYP19A1 active site were compared. Compounds 5b,c formed hydrophobic π−Sulfur interactions with the Hem600 structure. However, letrozole formed π−π T-shaped hydrophobic and πdonor hydrogen bond interactions with the Hem600 structure. Compounds 5a,b and and letrozole formed hydrophobic interactions with Phe221. Protein−ligand interaction details of letrozole with the CYP19A1 enzyme are given in Table 3.
2.5. Molecular Dynamics Simulations. Molecular dynamics (MD) simulations are frequently used to detect the in silico variation of protein−ligand interactions. The stability of the protein−ligand complex structure obtained from the molecular docking study can be estimated by calculating with an MD simulation. In this study, the protein− ligand complexes obtained from the molecular docking study of compounds 5b,c against the aromatase enzyme and the apoprotein without ligand were run under the same conditions for 100 ns to provide an MD validation. Topology files of the 5b,c ligands, the CYP19A1 enzyme, and the HEM structure it carries were obtained using the charmm36-feb2021 force field. RMSD and RMSF trajectory analyzes were performed. RMSD is a parameter that is frequently used to examine deviations and changes in protein structure. The RMSD measurement was measured by preferring the backbone atoms of the protein. CYP19A1 apoprotein, CYP19A1−5b, and CYP19A1−5c protein−ligand complex structures were measured as the means 0.32, 0.23, and 0.27 nm, respectively. As seen in Figure  2, the protein became more stable because of the interaction of 5b,c with CYP19A1. In particular, the CYP19A1−5b complex was stabilized below 0.2 nm. The other analysis parameter, RMSF, shows fluctuations and conformational changes in the protein's structure. A lower RMSF value indicates that the protein is more stable. As shown in Figure 2, residues at the Nand C-termini of the CYP19A1 apoprotein fluctuated more than the case for the CYP19A1−5b and CYP19A1−5c complex. Compounds 5b,c formed strong interactions with the aromatase enzyme, making CYP19A1 more stable. The 10000th frame protein−ligand interactions are presented in Figure 3 to explain the changes in the CYP19A1 active site at the end of 100 ns simulations of compounds 5b,c. The π−sulfur interaction of Hem600 with the thioether group of 5b,c was preserved. Compound 5b formed Asp309 with the benzimidazole core −NH and Gln218 with the phenol −OH. In compound 5c, H bonds were formed between the O atom of the oxazole ring and Thr310 and between fluorine and Trp224. As shown in Figure 3, 5b was extremely stable (mean RMSD: 0.14 nm) during the MD simulation, while 5c was stable up to about 40 nm, after which an increase to 0.4 nm (mean RMSD: 0.21 nm) was measured.
Another method to understand the stability of protein− ligand interactions in MD simulations is to measure protein− ligand interaction energies over time. The average short-range Coulombic interaction energy, the short-range Lennard−Jones energy, and MMPBSA for bindin -free energy (BFE) were measured for the protein−ligand interaction energy formed by compounds 5b,c by the CYP19A1 enzyme. BFE was obtained from the sum of van der Waals, electrostatic, polar solvation, and SASA energies. As seen in Table 4 compounds 5b,c formed strong interactions at the CYP19A1 enzyme active site with small deviations.
2.6. 3D-QSAR Study. In order to complement the experimental analysis, we carried out a 3D-QSAR study using a comparative molecular field analysis (CoMFA) and a comparative molecular similarity index analysis (CoMSIA). The study was performed according to the previous report of our research group. 18−22 Due to the scarce information on benzimidazole derivatives as aromatase inhibitors, we used the indole and benzofuran cores for the molecular alignment (all details are given in Table S1). All of the compounds span a biological activity range of more than 3 logarithmic units ( Figure S2). The best model of each analysis (CoMFA and CoMSIA) was selected by a combination of the following feature combination: q 2 > 0.5, a low number of components, and high value for the regression coefficient r 2 ncv . The models CoMFA-SE and CoMSIA-SH met these statistical criteria (Tables S2 and S3). A statistical summary of the best models (q 2 , N, SEP, SEE, r 2 ncv , F, PRESS, SD, and r 2 pred ), is reported in Table S4.
The aromatase inhibition values, with the predictions made by the selected CoMFA-SE and CoMSIA-SH models, are presented in Table S5. The experimental versus predictive activity graphs for the best CoMFA and CoMSIA models were created to visualize whether there is an adequate linear distribution of the predictive results for both models ( Figure  4). A good data distribution along the line y = x for both CoMFA and CoMSIA models for aromatase inhibitors in the training and test sets was observed. Moreover, each model was validated statistically using the Golbraikh and Tropsha methods. 23−26 The CoMFA and CoMSIA models reported here fulfilled the main conditions for internal and external validation (Table S6). Finally, to avoid a random correlation, a Y-random test was executed, giving a low value of q 2 (q 2 = −0.225 to 0.012) and r 2 ncv values below 0.5, confirming that both models developed are not a product of chance correlation (Table S7).
Finally, we performed the calculation of the applicability domain to corroborate the chemical suitability of all the compounds used in the studies. The applicability domain (AD) is a theoretical region in chemical space encompassing both the model descriptors and modeled response. It allows estimating the uncertainty in the prediction of a compound on the basis of how similar it is to the training compounds employed in the model. In this work, we used the method developed by Roy et al. for the determination of AD which is based on the basic theory of the standardization approach. 22 In our cases, all the training and test set molecules, from both CoMFA and CoMSIA, were found within the applicability domain.
2.6.1. Contour Map Analysis. The advantage of the use of the 3D-QSAR technique is obtaining contour maps around the target molecule. These contour maps show the main features as being favorable or unfavorable for the biological activity on the basis of the different colored polyhedra around a molecule. To visualize the contour maps, the compound 2-((1Himidazol-1-yl)methyl)-1-(4-(trifluoromethyl)phenyl)-1H-indole (13) was used in the CoMFA-SE and CoMSIA-SH diagrams ( Figure 5).
The CoMFA and CoMSIA steric maps of aromatase inhibitors are very similar ( Figure 5A−C). In both figures, we can see a green polyhedron around position 6 of the indole core. This means that the bulk substituent favors the inhibitory activity. In fact, several indole derivatives with bulky substituents on position 6 have pIC 50 > 6.0, for example, the compounds 11−28 (Me, OMe, or Cl substituents), as well as a similar substituent on position 6 of benzofuran core as in compounds 89−99. Moreover, the substituent on position 5 of the indole core can project some fragments to near the green polyhedron, as for compound 63 with the phenyl group bonded to amino-sulfonyl. Interestingly, compounds 1−8 have a bulky substituent on chiral carbon that can project the Nimidazole or phenyl rings toward the green polyhedra. This phenomenon favors aromatase inhibition activity, concordant with pIC 50 > 6.0. Moreover, a bulky substitution on the Natom of the indole core improves the aromatase inhibition activity. This is due to the substituent being near the green polyhedron ( Figure Figure  5B) shows two main areas: the first between position 7 and the N-substituent of the indole core (red polyhedron) and the second between the N-substituent and position 2 of the indole core (blue polyhedron). In the first case, an electron-rich substituent improves the aromatase inhibitory activity as in compounds 108, 109, and 112−115 (pIC 50 > 6.0). All of them have a halogen atom oriented toward the red polyhedron. In the second case, an electron-rich substituent causes a decrease in the aromatase inhibitory activity (pIC 50 < 6.0), as in compounds 47−51 and 53−62. These compounds have a 2indolin substituent bonded to position 2 of the indole core, oriented to the blue polyhedron.
The hydrophobic map of the CoMSIA model ( Figure 5D) shows a white polyhedron near positions 5 and 6 of the indole core; this means that a hydrophobic substituent decreases the aromatase activity as in compounds 9, 66, 68, 70, 72, 79−81, 100, 122, and 123, which have a methyl, halogen, or methoxy group and other substituents, as well as a value of pIC 50 < 6.0. A similar effect is shown between positions 2 and 3 of the indole core, where a hydrophobic substituent causes a decrease in aromatase inhibitory activity as in compounds 29−45, 56− 62, 64−66, 68, 70, 72, 100−110, 113−116, 122, and 123, which have analkyl or aromatic ring bonded to an alkyl group. Moreover, between the N-substituent and the group bonded to position 2 of the indole core, there is a yellow polyhedron, which means that the hydrophobic group causes an increased aromatase inhibition activity (pIC 50 > 6.0). Some of the compounds with a hydrophobic group in this position are 11− 28, 67, 69, 71, 73−78, 89−99, 108−112, and 117−121 (e.g., N-alkyl, N-aryl, or an alkyl-aromatic substituent bonded to position 2 of indole core).
2.6.2. Experimental Validation. To externally validate the models, we carried out the synthesis of four benzimidazoletype molecules. These compounds were predicted to have biological activity by the models reported here, and their experimental activity was determined by enzymatic assays. The results obtained are shown in Table 5.
The predictive values based on CoMFA and CoMSIA models of aromatase inhibitors showed higher values in comparison to the experimental values. To explain this phenomenon, we overlaid compound 5b on the contour map obtained from the CoMFA and CoMSIA study ( Figure 6). In the steric map, the phenyl ring bonded to the benzimidazole moiety is oriented toward the yellow polyhedra, causing a decrease in the aromatase inhibition activity ( Figure 6A−C). Moreover, the 1,3,4-oxadiazole fragment is projected near the green polyhedra, favoring the inhibition activity on aromatase ( Figure 6A−C). In the case of the electrostatic map, no substituents on 5b are close to the red or blue polyhedron ( Figure 6B), and so there is no increase or decrease in aromatase inhibition activity. Finally, in the hydrophobic map ( Figure 6D), the 1,3,4-oxadiazole fragment bonded to the benzimidazole core is close to the yellow polyhedra, favoring the aromatase inhibition activity; however, the phenyl ring bonded to position 2 of benzimidazole core is mainly near the white polyhedron, causing a decrease in\aromatase inhibition activity.

CONCLUSION
In the present investigation, benzimidazole-1,3,4-oxadiazole derivatives 5a−i have been synthesized and evaluated for their in vitro anticancer properties. The structures of synthesized compounds 5a−i were characterized by using various modern analytical techniques such as 1 H NMR, 13 C NMR, and HRMS. Among the compounds, 5b,c were found to be the most potent compounds against the MCF-7 cell line. According to the result of in vitro aromatase inhibitory activity, compounds 5b,c cause 50% aromatase enzyme inhibition effectively at 1.475 ± 0.062 and 2.512 ± 0.124 μM, respectively. Molecular docking studies were performed to find the interaction types with the aromatase enzyme for compounds 5b,c. Furthermore, a molecular dynamics simulation was performed out to explore the most likely binding mode of compounds 5b,c with CYP19A1. Two types of satisfactory 3D-QSAR (CoMFA and CoMSIA) models were generated, to predict the inhibitory activities of novel inhibitors. In our cases, all of the training and test set molecules, from both CoMFA and CoMSIA, were found to be within the applicability domain.  In DMF, a combination of a 4-substituted benzaldehyde (0.03 mol) and sodium metabisulfide (0.03 mol, 5.7 g) was irradiated with microwaves for 5 min at 240°C and 10 bar (Anton-Paar Monowave 300). After the mixture had cooled, 3,4-diaminobenzoic acid (0.03 mol, 4.56 g) was added and the reaction conditions were maintained. TLC was used to track the reaction's progress. The mixture was put into crushed ice after the reaction was completed.

Synthesis of Methyl 2-(4-Substituted phenyl)-1Hbenzo[d]
imidazole-6-carboxylate derivatives 2a−c. In methanol, were dissolved compounds 1a−c (0.025 mol), and several drops of sulfuric acid were added. The mixture was refluxed for 72 h before filtering the precipitate. 16 4.2.4. Synthesis of 2-(4-Substituted phenyl)-1H-benzo[d]imidazole-6-carbohydrazide derivatives 3a−c. In a microwave synthesis reactor vial (30 mL), compounds 2a−c (0.018 mol) were combined with an excess of hydrazine hydrate (5 mL) in EtOH (15 mL) and irradiated with microwaves (Anton-Paar Monowave 300). The reaction mixture was heated for 10 min at 240°C and 10 bar. The mixture was cooled and then poured into ice water, filtered, rinsed with water, dried, and recrystallized from EtOH. 16 4.2.5. Synthesis of 2-((4-Substituted phenyl)-(6-(5-mercapto-1,3,4-oxadiazol-2-yl)-1H-benzo[d]imidazole Derivatives 4a−c. In a solution of NaOH (0.01 mol, 0.4 g) in  ethanol were dissolved compounds 3a−c (0.01 mol), and carbon disulfide (0.01 mol, 0.60 mL) was added to the mixture. The mixture was refluxed for 8 h. The solution was cooled and acidified to pH 4−5 using concentrated hydrochloric acid solution after the reaction was completed to obtain the target compound. 16 4.2.6. General Synthesis Method of Target Compounds 5a−i. In acetone were dissolved compounds 4a−c (0.001 mol), potassium carbonate (0.001 mol, 0.138 g), and 2-chloro-1-(4-substituted epiperazin-1-yl)ethane-1-one (0.0015 mol) and refluxed for 6 h. The solvent was evaporated after TLC control, the residue was washed with water, dried, and then recrystallized from ethanol to produce final compounds 5a− i. 16 4.2.6.1. 2-((5-(2-(4-Hydroxyphenyl)-1H-benzo[d]imidazol-6-yl)-1,3,4-oxadiazol-2-yl)thio)-1-(4-benzylpiperazin-1-yl)ethan-1-one (5a  29 The topology file of CYP19A1 including the HEM structure was created with the Charmm36 force field, and ligand topology files were created via the CHARMM General Force Field (CGenFF) server. 30 The topology files for CYP19A1 and compounds 5b,c were combined and solved with the TIP3 water model, and the appropriate Na + and Cl − ions were added. The energy of the system that formed CYP19A1, HEM, ligand, ion, and solvent was minimized. The system was equilibrated with 0.2 ns NVT and 0.3 ns NPT stages at 1 atm pressure and 300 K temperature according to a V-rescale thermostat and Berendsen barostat, respectively, and 100 ns molecular dynamics simulations were performed to 2 fs. The trajectory analysis was performed with gmx scripts, and the root-mean-square deviation (RMSD) and the root-meansquare fluctuation (RMSF) were measured. The trajectory results were visualized with VMD (Visual Molecular Dynamics), and RMSD and RMSF graphs were created with GraphPad Prism 8 software. The binding free energy calculation by the molecular mechanics Poisson−Boltzmann surface area (MM-PBSA) was calculated between 80 and 100 ns using RashmiKumari's g_mmpbsa package. 31 4.6. 3D-QSAR Study. 3D-QSAR studies were performed with Sybyl X-1.2 software installed in a Windows 10 environment on a PC with an Intel core i7 CPU, according to previous reports of our research group. 18−20 The details of the procedure, data set, statistical values, and the validation of models are given in the Supporting Information.

* sı Supporting Information
The Supporting Information is available free of charge at It includes the The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acsomega.2c01497. 1 H NMR, 13